function main()

[data,labels]=loadDigits(2000,'train');
j_0=1;
j_1=1;
j_2=1;
j_3=1;
j_4=1;
j_5=1;
j_6=1;
j_7=1;
j_8=1;
j_9=1;
for i=1:2000
    if labels(1,i)==0
        a0(:,j_0)=data(:,i);
        j_0=j_0+1;
    end
    if labels(1,i)==1
        a1(:,j_1)=data(:,i);
        j_1=j_1+1;
    end
    if labels(1,i)==2
        a2(:,j_2)=data(:,i);
        j_2=j_2+1;
    end
    if labels(1,i)==3
        a3(:,j_3)=data(:,i);
        j_3=j_3+1;
    end
    if labels(1,i)==4
        a4(:,j_4)=data(:,i);
        j_4=j_4+1;
    end
    if labels(1,i)==5
        a5(:,j_5)=data(:,i);
        j_5=j_5+1;
    end
    if labels(1,i)==6
        a6(:,j_6)=data(:,i);
        j_6=j_6+1;
    end
    if labels(1,i)==7
        a7(:,j_7)=data(:,i);
        j_7=j_7+1;
    end
    if labels(1,i)==8
        a8(:,j_8)=data(:,i);
        j_8=j_8+1;
    end
    if labels(1,i)==9
        a9(:,j_9)=data(:,i);
        j_9=j_9+1;
    end
end
[m1,n1]=size(a1);
[m4,n4]=size(a4);
[m7,n7]=size(a7);
[m9,n9]=size(a9);
for i=1:n1
    a17(:,i)=a1(:,i);
    label17(i)=1;
end
for i=n1+1 : n1+n7
    a17(:,i)=a7(:,i-n1);
    label17(i)=7;
end
for i=1:n4
    a49(:,i)=a4(:,i);
    label49(i)=4;
end
for i=n4+1 : n4+n9
    a49(:,i)=a9(:,i-n4);
    label49(i)=9;
end
d=tangent_d(data,data);
d17=tangent_d(a17,a17);
d49=tangent_d(a49,a49);
options.comp=1;
options.display=1;
options.verbose=1;
options.overlay=1;
options.dims=1:10;
[Y,R,E]=Isomap(d,'k',10,options);
[Y17,R17,E17]=Isomap(d17,'k',10,options);
[Y49,R49,E49]=Isomap(d49,'k',10,options);

dims=options.dims;
comp = options.comp; 
overlay = options.overlay; 
displ = options.display; 
verbose = options.verbose;

 figure;
     hold on
     plot(dims, R, 'bo'); 
     plot(dims, R, 'b-'); 
     hold off
     ylabel('Residual variance'); 
     xlabel('Isomap dimensionality'); 

     %%%%% Plot two-dimensional configuration %%%%%
     twod = find(dims==2); 
     if ~isempty(twod)
         figure;
         hold on;
         for i=1:2000
             if labels(1,i)==1 
                plot(Y.coords{twod}(1,i), Y.coords{twod}(2,i), 'ro'); 
             end
             if labels(1,i)==2
                 plot(Y.coords{twod}(1,i), Y.coords{twod}(2,i), 'bo');
             end
             if labels(1,i)==3
                 plot(Y.coords{twod}(1,i), Y.coords{twod}(2,i), 'yo');
             end
             if labels(1,i)==4
                 plot(Y.coords{twod}(1,i), Y.coords{twod}(2,i), 'mo');
             end
             if labels(1,i)==5
                 plot(Y.coords{twod}(1,i), Y.coords{twod}(2,i), 'co');
             end
             if labels(1,i)==6
                 plot(Y.coords{twod}(1,i), Y.coords{twod}(2,i), 'go');
             end
             if labels(1,i)==7
                 plot(Y.coords{twod}(1,i), Y.coords{twod}(2,i), 'ko');
             end
             if labels(1,i)==8
                 plot(Y.coords{twod}(1,i), Y.coords{twod}(2,i),'m*');
             end
             if labels(1,i)==9
                 plot(Y.coords{twod}(1,i), Y.coords{twod}(2,i),'rx');
             end
             if labels(1,i)==0
                 plot(Y.coords{twod}(1,i), Y.coords{twod}(2,i), 'b+');
             end
         end
         legend('0','1','2','3','4','5','6','7','8','9');
         if (overlay == 1)
             %gplot(E(Y.index, Y.index), [Y.coords{twod}(1,:); Y.coords{twod}(2,:)]'); 
             title('Two-dimensional Isomap embedding (with neighborhood graph).'); 
         else
             title('Two-dimensional Isomap.'); 
         end
         hold off;
     end
     
     Imap(Y17,R17,1,7,options,label17,n1+n7,a17);
     Imap(Y49,R49,4,9,options,label49,n9+n4,a49);
